datos_tp <-read_excel("DatosTP1.xlsx") %>% mutate_at('variedad', as.factor)
glimpse(datos_tp)
## Rows: 6,497
## Columns: 13
## $ `acidez fija` <dbl> 7.0, 6.3, 8.1, 7.2, 7.2, 8.1, 6.2, 7.0, 6.…
## $ `acidez volátil` <dbl> 0.27, 0.30, 0.28, 0.23, 0.23, 0.28, 0.32, …
## $ `ácido cítrico` <dbl> 0.36, 0.34, 0.40, 0.32, 0.32, 0.40, 0.16, …
## $ `azúcar residual` <dbl> 20.70, 1.60, 6.90, 8.50, 8.50, 6.90, 7.00,…
## $ cloruros <dbl> 0.045, 0.049, 0.050, 0.058, 0.058, 0.050, …
## $ `anhídrido sulfuroso libre` <dbl> 45, 14, 30, 47, 47, 30, 30, 45, 14, 28, 11…
## $ `anhídrido sulfuroso total` <dbl> 170, 132, 97, 186, 186, 97, 136, 170, 132,…
## $ densidad <dbl> 1.0010, 0.9940, 0.9951, 0.9956, 0.9956, 0.…
## $ pH <dbl> 3.00, 3.30, 3.26, 3.19, 3.19, 3.26, 3.18, …
## $ sulfatos <dbl> 0.45, 0.49, 0.44, 0.40, 0.40, 0.44, 0.47, …
## $ alcohol <dbl> 8.8, 9.5, 10.1, 9.9, 9.9, 10.1, 9.6, 8.8, …
## $ calidad <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5, 5, 7, …
## $ variedad <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
set.seed(661) #fijo la semilla
sample_var_1 = datos_tp %>% filter(datos_tp$variedad == 1)
sample_var_2 = datos_tp %>% filter(datos_tp$variedad == 2)
index_sample_data_1 <- sample(1:nrow(sample_var_1), size=1000, replace=FALSE)
index_sample_data_2 <- sample(1:nrow(sample_var_2), size=1000, replace=FALSE)
sample_data_var_1 = sample_var_1[index_sample_data_1,]
sample_data_var_2 = sample_var_2[index_sample_data_2,]
sample_data <- rbindlist(list(sample_data_var_1, sample_data_var_2)) # Rbind data.tables
glimpse(sample_data)
## Rows: 2,000
## Columns: 13
## $ `acidez fija` <dbl> 5.9, 5.8, 9.1, 6.9, 7.2, 7.2, 7.4, 5.2, 7.…
## $ `acidez volátil` <dbl> 0.320, 0.300, 0.280, 0.250, 0.230, 0.230, …
## $ `ácido cítrico` <dbl> 0.26, 0.27, 0.49, 0.35, 0.32, 0.19, 0.27, …
## $ `azúcar residual` <dbl> 1.5, 1.7, 2.0, 9.2, 8.5, 13.7, 1.2, 1.1, 1…
## $ cloruros <dbl> 0.057, 0.014, 0.059, 0.034, 0.058, 0.052, …
## $ `anhídrido sulfuroso libre` <dbl> 17.0, 45.0, 10.0, 42.0, 47.0, 47.0, 27.0, …
## $ `anhídrido sulfuroso total` <dbl> 141, 104, 112, 150, 186, 197, 99, 69, 112,…
## $ densidad <dbl> 0.99170, 0.98914, 0.99580, 0.99470, 0.9956…
## $ pH <dbl> 3.24, 3.40, 3.15, 3.21, 3.19, 3.12, 3.19, …
## $ sulfatos <dbl> 0.36, 0.56, 0.46, 0.36, 0.40, 0.53, 0.33, …
## $ alcohol <dbl> 10.700000, 12.600000, 10.100000, 11.500000…
## $ calidad <dbl> 5, 7, 5, 6, 6, 5, 6, 6, 5, 5, 7, 6, 6, 6, …
## $ variedad <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
print(paste("Balanced dataset - Variedad 1: ", nrow((sample_data %>% filter(variedad == 1)))))
## [1] "Balanced dataset - Variedad 1: 1000"
print(paste("Balanced dataset - Variedad 2: ", nrow((sample_data %>% filter(variedad == 2)))))
## [1] "Balanced dataset - Variedad 2: 1000"
colnames(sample_data) <- c('acidez_fija','acidez_volatil','acido_citrico','azucar_residual','cloruros','anhidrido_sulfuroso_libre','anhidrido_sulfuroso_total','densidad','pH','sulfatos','alcohol','calidad','variedad')
sample_data
## acidez_fija acidez_volatil acido_citrico azucar_residual cloruros
## 1: 5.9 0.32 0.26 1.5 0.057
## 2: 5.8 0.30 0.27 1.7 0.014
## 3: 9.1 0.28 0.49 2.0 0.059
## 4: 6.9 0.25 0.35 9.2 0.034
## 5: 7.2 0.23 0.32 8.5 0.058
## ---
## 1996: 5.6 0.54 0.04 1.7 0.049
## 1997: 7.3 0.65 0.00 1.2 0.065
## 1998: 6.9 0.41 0.31 2.0 0.079
## 1999: 9.0 0.38 0.41 2.4 0.103
## 2000: 9.2 0.67 0.10 3.0 0.091
## anhidrido_sulfuroso_libre anhidrido_sulfuroso_total densidad pH
## 1: 17 141 0.99170 3.24
## 2: 45 104 0.98914 3.40
## 3: 10 112 0.99580 3.15
## 4: 42 150 0.99470 3.21
## 5: 47 186 0.99560 3.19
## ---
## 1996: 5 13 0.99420 3.72
## 1997: 15 21 0.99460 3.39
## 1998: 21 51 0.99668 3.47
## 1999: 6 10 0.99604 3.13
## 2000: 12 48 0.99888 3.31
## sulfatos alcohol calidad variedad
## 1: 0.36 10.7 5 1
## 2: 0.56 12.6 7 1
## 3: 0.46 10.1 5 1
## 4: 0.36 11.5 6 1
## 5: 0.40 9.9 6 1
## ---
## 1996: 0.58 11.4 5 2
## 1997: 0.47 10.0 7 2
## 1998: 0.55 9.5 6 2
## 1999: 0.58 11.9 7 2
## 2000: 0.54 9.5 6 2
df_pca = sample_data
df_numericas_pca <-df_pca %>% dplyr::select(where(is.numeric))
data_long <- melt(df_pca)
## Using variedad as id variables
#data_long
ggplot(data_long, aes(x=variable, y=value)) +
geom_boxplot() +
facet_wrap(~variable, scale="free")
#Boxplot diferenciado por variedad de vino
ggplot(data_long, aes(x=variable, y=value, fill= variedad)) +
geom_boxplot() +
facet_wrap(~variable, scale="free")
#Matriz de correlación
m_cor <- cor(df_numericas_pca)
# representa la matriz de correlaciones mediante círculos
corrplot(m_cor,method="circle")
variance_pca <- df_pca %>%
summarise_if(is.numeric, var) %>%
t() %>%
data.frame() %>%
round(3)
variance_pca
## .
## acidez_fija 2.436
## acidez_volatil 0.036
## acido_citrico 0.028
## azucar_residual 17.975
## cloruros 0.002
## anhidrido_sulfuroso_libre 339.797
## anhidrido_sulfuroso_total 3715.930
## densidad 0.000
## pH 0.027
## sulfatos 0.030
## alcohol 1.352
## calidad 0.761
pca <- prcomp(df_numericas_pca,scale = TRUE)# con datos estandarizados
names(pca)
## [1] "sdev" "rotation" "center" "scale" "x"
round(pca$rotation,2)
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9
## acidez_fija 0.27 -0.31 0.38 -0.25 0.14 -0.17 0.35 -0.20 0.20
## acidez_volatil 0.39 0.02 -0.31 -0.03 0.06 0.20 0.61 0.10 0.10
## acido_citrico -0.13 -0.27 0.51 -0.02 -0.15 -0.40 -0.06 0.38 0.13
## azucar_residual -0.31 -0.34 -0.14 0.09 0.47 0.19 0.03 0.43 -0.38
## cloruros 0.31 -0.24 0.02 0.36 -0.42 0.40 -0.14 0.45 0.27
## anhidrido_sulfuroso_libre -0.40 -0.13 -0.13 0.38 -0.16 -0.13 0.42 -0.18 0.33
## anhidrido_sulfuroso_total -0.46 -0.15 -0.12 0.19 -0.16 -0.04 0.22 -0.05 0.01
## densidad 0.22 -0.51 -0.11 0.07 0.42 -0.13 0.00 0.00 0.06
## pH 0.20 0.30 -0.31 0.35 0.26 -0.59 -0.17 0.22 0.27
## sulfatos 0.29 -0.10 0.22 0.63 -0.08 -0.10 0.03 -0.35 -0.55
## alcohol -0.01 0.46 0.33 0.06 0.06 -0.06 0.47 0.43 -0.25
## calidad -0.12 0.22 0.41 0.31 0.50 0.43 -0.09 -0.16 0.41
## PC10 PC11 PC12
## acidez_fija -0.29 -0.32 -0.43
## acidez_volatil 0.52 0.17 -0.07
## acido_citrico 0.47 0.27 0.02
## azucar_residual -0.09 0.09 -0.41
## cloruros -0.22 -0.16 -0.04
## anhidrido_sulfuroso_libre -0.34 0.44 -0.01
## anhidrido_sulfuroso_total 0.33 -0.72 0.08
## densidad -0.10 -0.05 0.68
## pH -0.01 -0.16 -0.24
## sulfatos 0.11 0.05 -0.09
## alcohol -0.29 -0.13 0.32
## calidad 0.20 -0.02 0.00
contrib <- as.matrix(round(pca$rotation,2))
corrplot(contrib,is.corr=FALSE)
round(pca$center,2) #vector de medias
## acidez_fija acidez_volatil acido_citrico
## 7.60 0.40 0.31
## azucar_residual cloruros anhidrido_sulfuroso_libre
## 4.52 0.07 25.93
## anhidrido_sulfuroso_total densidad pH
## 93.24 1.00 3.25
## sulfatos alcohol calidad
## 0.58 10.47 5.76
prop_varianza <- pca$sdev^2 / sum(pca$sdev^2)
prop_varianza #varianza explicada por cada componente
## [1] 0.27683825 0.20521757 0.15634119 0.08242453 0.06948051 0.04950167
## [7] 0.04356408 0.03971796 0.03753468 0.02186397 0.01427615 0.00323943
prop_varianza_acum_pca <- cumsum(prop_varianza)
prop_varianza_acum_pca #varianza acumilada explicada por los componentes de forma creciente segun su importancia
## [1] 0.2768383 0.4820558 0.6383970 0.7208215 0.7903021 0.8398037 0.8833678
## [8] 0.9230858 0.9606205 0.9824844 0.9967606 1.0000000
ggplot(data = data.frame(prop_varianza_acum_pca, pc = 1:12),
aes(x = pc, y = prop_varianza_acum_pca, group = 1)) +
geom_point() +
geom_line() +
theme_bw() +
labs(x = "Componente principal",
y = "Prop. varianza explicada acumulada")
screeplot(pca, type = "l", npcs = 12)
abline(h = 1, col="red", lty=5)
legend("topright", legend=c("Eigenvalue = 1"),
col=c("red"), lty=5, cex=0.6)
fviz_eig(pca, ncp =12, addlabels = TRUE, main="")
#pca_2$x[,1]# scores para la primer componente (PC1)
res.ind <- get_pca_ind(pca)
head(res.ind$coord,10)
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6
## 1 -1.01189442 1.3191155 -0.77274712 -1.0357319 -1.303163205 -0.04051432
## 2 -1.66875699 3.3360841 0.67248805 0.9158157 0.003178629 -0.32182485
## 3 0.01602419 -0.8192638 0.82974834 -1.5371725 -0.683656840 -0.75875675
## 4 -2.33345523 0.2698171 0.04272720 -0.2733312 0.534078552 -0.07713830
## 5 -2.33834010 -0.7957860 -0.42077524 0.1231915 0.134039028 0.15845606
## 6 -2.23501071 -2.3477902 -1.56768877 0.2210427 0.509138131 0.22625340
## 7 -1.14291956 0.7835540 0.05342973 -1.2018491 -0.456717256 0.13685204
## 8 -0.71628256 1.6486190 -0.14371084 -0.3971398 -0.535746151 0.10267295
## 9 -0.51483343 0.4081233 -0.23681707 -1.5712144 -1.066796279 0.21658374
## 10 -0.68431517 1.1294388 0.37372431 -0.7401340 -0.808225877 -1.52074829
## Dim.7 Dim.8 Dim.9 Dim.10 Dim.11 Dim.12
## 1 -0.495919296 0.35113630 -0.05141232 0.25368656 -0.6094454 0.150398260
## 2 0.470909358 -0.13683111 0.23585111 -0.06141238 0.3736711 -0.280061337
## 3 -0.391449922 0.06641265 0.11350261 0.21009852 -0.6343479 0.119948769
## 4 0.435569228 0.76019151 0.04545350 -0.41840389 -0.1270732 0.180397465
## 5 -0.009342246 0.04604676 0.55799854 -0.22347452 -0.4765575 -0.015741529
## 6 -0.035025845 -0.30334927 -0.77677646 -0.63039045 -0.5180117 0.009249802
## 7 -0.815847371 -0.65442870 0.89109712 -0.22575357 -0.1096548 -0.145302705
## 8 -1.296704575 -0.60975149 -0.23873680 0.62686972 0.6294417 0.047023266
## 9 -0.786231545 -0.47653397 -0.34558505 -0.01493300 -0.5615664 -0.048033200
## 10 -1.155540815 0.44736925 -0.28989256 0.03669330 -0.1287600 -0.015182979
autoplot(pca,
data = df_pca,
colour = 'variedad',
loadings = TRUE,
loadings.colour = 'black',
loadings.label = TRUE,
loadings.label.size = 3)
Conclusiones
Podemos obtener como resultado del PCA, que sería prudente quedarnos con las primeras 7 componentes que representan más del 86% de la variabilidad de la información, además se comienza a aplanar la curva que explica la diferencia de vairabilidad entre las siguientes componentes
Observamos que los principales atributos con los que nos quedamos a posteriri de decidir con la información de las componentes serian anhidrido_sulfuroso_total (1) densidad (2) acido_citrico (3) sulfatos (4) calidad (5) pH (6) alcohol (7)
anhidrido_sulfuroso_total y anhidrido_sulfuroso_libre correlacionan de forma casi total
Mientras que otras variables tienne una correlación alta dada por el ángulo que forman sus autovectores Azucar_resicual/Acido_citrico, acidez_fija/cloruros, sulfatos/acidez_volatil
Por otro lado, vemos que otras variables no tienen ningún tipo de correlación y las podemos determinar independientes entre si, como alcohol/acidez_volatil, ph/calidad, acidez_fija/azucar_residual
#Separo el dataset para trabajar el analisis discriminante
df_ad = sample_data
ggpairs(df_ad, legend = 1, columns = 1:12, aes(color = variedad, alpha = 0.5),
upper = list(continuous = "points"))+
theme(legend.position = "bottom")
#Separo en train y test para predecir
set.seed(661)#setear la semilla
# Create data split for train and test
df_split_ad <- initial_split(df_ad,
prop = 0.9,
strata = variedad)#en este caso no es necesario porque las clases están balanceadas
df_train_ad <- df_split_ad %>% training()
df_test_ad <- df_split_ad %>% testing()
# Número de datos en test y train
paste0("Total del dataset de entrenamiento: ", nrow(df_train_ad))
## [1] "Total del dataset de entrenamiento: 1800"
paste0("Total del dataset de testeo: ", nrow(df_test_ad))
## [1] "Total del dataset de testeo: 200"
blanco_ad<- subset(df_ad[,1:12], df_ad$variedad == "1")
tinto_ad <- subset(df_ad[,1:12], df_ad$variedad == "2")
mvShapiro.Test(as.matrix(blanco_ad))
##
## Generalized Shapiro-Wilk test for Multivariate Normality by
## Villasenor-Alva and Gonzalez-Estrada
##
## data: as.matrix(blanco_ad)
## MVW = 0.91982, p-value < 2.2e-16
#Valido normalidad => no se cumple porque es < 0,05
mvShapiro.Test(as.matrix(tinto_ad))
##
## Generalized Shapiro-Wilk test for Multivariate Normality by
## Villasenor-Alva and Gonzalez-Estrada
##
## data: as.matrix(tinto_ad)
## MVW = 0.88076, p-value < 2.2e-16
#Matriz vairanza y covarianza son iguales, valido homocedasticidad => no se cumple tampoco
boxM(data = df_ad[, 1:12], grouping = df_ad$variedad)
##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: df_ad[, 1:12]
## Chi-Sq (approx.) = 4289.5, df = 78, p-value < 2.2e-16
#Por mas que nada se cumple seguimos el analisis, sabiendo que los resultados pueden no tener la validez que esperamos
Conclusiones
Determino independencia suponiendo que los datos se obtuvieron de fuentes independientes
Vemos que no se cumple la normalidad en ninguna de las dos distribuciones planteadas
Vemos que no se cumple homocedasticidad en el conjunto de datos
Por mas que los supuestos no se cumplen, avanzamos con el análisis a sabiendas que los resultados pueden no tener la validez que esperamos
model_lda <- lda(variedad~., data=df_train_ad)
model_lda
## Call:
## lda(variedad ~ ., data = df_train_ad)
##
## Prior probabilities of groups:
## 1 2
## 0.5 0.5
##
## Group means:
## acidez_fija acidez_volatil acido_citrico azucar_residual cloruros
## 1 6.860278 0.2821444 0.3411222 6.524278 0.04621778
## 2 8.372889 0.5232333 0.2759889 2.574833 0.08808333
## anhidrido_sulfuroso_libre anhidrido_sulfuroso_total densidad pH
## 1 36.08000 140.10167 0.9941380 3.186922
## 2 15.94556 46.83222 0.9967918 3.305744
## sulfatos alcohol calidad
## 1 0.4923333 10.47361 5.858889
## 2 0.6612444 10.43659 5.646667
##
## Coefficients of linear discriminants:
## LD1
## acidez_fija -0.44594457
## acidez_volatil 2.05824631
## acido_citrico -0.58935312
## azucar_residual -0.36366284
## cloruros 2.40213219
## anhidrido_sulfuroso_libre 0.01517737
## anhidrido_sulfuroso_total -0.01876167
## densidad 898.41125911
## pH -1.19965511
## sulfatos 0.92102307
## alcohol 0.75018119
## calidad 0.05713148
prop = model_lda$svd^2/sum(model_lda$svd^2)
prop
## [1] 1
#Coeficiente lda calculado para cada entrada del set analizado
lda.data <- cbind(df_train_ad, predict(model_lda)$x)
#lda.data
#ggplot(lda.data, aes(LD1)) +
# geom_point(aes(color = variedad))
#Predicciones sobre el set de entrenamiento para cada entrada del set analizado
predictions_train <- model_lda %>% predict(df_train_ad)
#predictions_train$class
table(predict(model_lda,type="class")$class,df_train_ad$variedad)
##
## 1 2
## 1 898 12
## 2 2 888
#partimat (variedad~. , data=df_train_ad , method="lda")
#ggord(model_lda, df_train_ad$variedad)
predictions <- model_lda %>% predict(df_test_ad)
predictions$class
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2
## [112] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [149] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [186] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## Levels: 1 2
lda.test <- predict(model_lda,df_test_ad)
df_test_ad$lda <- lda.test$class
table(df_test_ad$lda,df_test_ad$variedad)
##
## 1 2
## 1 100 0
## 2 0 100
Si los supuesto no se cumplen, podrían no ser válidos los resultados obtenidos Como se utiliza las matrices de varianzas y covarianzas, estas técnicas son sensibles a outliers
model_qda <- qda(variedad ~ ., df_train_ad)
model_qda
## Call:
## qda(variedad ~ ., data = df_train_ad)
##
## Prior probabilities of groups:
## 1 2
## 0.5 0.5
##
## Group means:
## acidez_fija acidez_volatil acido_citrico azucar_residual cloruros
## 1 6.860278 0.2821444 0.3411222 6.524278 0.04621778
## 2 8.372889 0.5232333 0.2759889 2.574833 0.08808333
## anhidrido_sulfuroso_libre anhidrido_sulfuroso_total densidad pH
## 1 36.08000 140.10167 0.9941380 3.186922
## 2 15.94556 46.83222 0.9967918 3.305744
## sulfatos alcohol calidad
## 1 0.4923333 10.47361 5.858889
## 2 0.6612444 10.43659 5.646667
table(predict(model_qda,type="class")$class,df_train_ad$variedad)
##
## 1 2
## 1 882 12
## 2 18 888
lda.test_qda <- predict(model_qda,df_test_ad)
df_test_ad$qda <- lda.test_qda$class
table(df_test_ad$qda,df_test_ad$variedad)
##
## 1 2
## 1 100 0
## 2 0 100
Conclusiones Podemos determinar a partir de la predicción y el analisis de su matriz de confusión, la separación generada por LDA en el conjunto de datos, arroja un accuracy del 100% en el dataset de prueba
Solo tenemos una dimensión debido a que nuestras variables objetivos solamente tienen dos posibilidades
Por el tipo de dato a validar y el dataset analizado, podría tener mayor sentido analizar otra variable para separa el mismo como por ejemplo la calidad donde nos arrojarían varias dimensiones a separar
En principio para el objetivo del análisis, determinamos que LDA funciona muy bien en su predicción de vinos tintos y blancos
No vemos una mejora significativa utilizando QDA
df_svm = sample_data
#Veo como correlacionan los datos
df_svm %>% ggpairs(.,mapping=ggplot2::aes(color = variedad,alpha = 0.1),
upper = list(continuous = wrap("cor", size = 2.5),discrete = "blank", combo="blank"),
lower = list(combo = "box"),progress = F)+
theme+
labs(title= 'Descripción de variables en la base de datos', x='Variable', y='Variable')+
scale_fill_manual(values=c('royalblue2','red'))+
scale_color_manual(values=c('royalblue2','#ff7474ff'))
# Grafico biplot de pca para poder comparar posteriormente con analisis de svm
datos_para_acp = df_svm %>%
dplyr::select(is.numeric) # todas las variables numéricas
## Warning: Predicate functions must be wrapped in `where()`.
##
## # Bad
## data %>% select(is.numeric)
##
## # Good
## data %>% select(where(is.numeric))
##
## ℹ Please update your code.
## This message is displayed once per session.
df_svm.pc = prcomp(datos_para_acp,scale = TRUE) #escalo los datos
#grafico
ggbiplot(df_svm.pc, obs.scale=0.1 ,var.scale=1,alpha=0.5,
groups=factor(df_svm$variedad)) +
scale_color_manual(name="variedad",values=c('royalblue2','#ff7474ff'),
labels=c("blanco",'tinto')) +
theme + labs(title='Análisis de componentes principales') +
theme(legend.position=c(.85,.15))
#Analisis estadistico univariado agrupado por variedad => Estaria bueno tenerlo pero no lo puedo hacer andar de momento
# Creo dataframes de acuerdo al diagnóstico
data_n <- df_svm%>%filter(variedad==1)
data_m <- df_svm%>%filter(variedad==2)
pval_all <- mshapiro.test(t(df_svm[,1:12]))
pval_blanco <- mshapiro.test(t(df_svm%>%dplyr::filter(variedad==1)%>% dplyr::select(1:12)))
pval_tinto <- mshapiro.test(t(df_svm%>%dplyr::filter(variedad==2)%>% dplyr::select(1:12)))
ShapiroW_p.valor <- c(pval_all$p.value, pval_blanco$p.value, pval_tinto$p.value)
Datos <- c('todos','subset blanco', 'subset tinto')
# Imprimo resultados de normalidad multivariada total y por grupos (estos últimos son relevantes)
kable(cbind(Datos,ShapiroW_p.valor))
| Datos | ShapiroW_p.valor |
|---|---|
| todos | 5.24820022494157e-43 |
| subset blanco | 1.32661116153684e-33 |
| subset tinto | 6.99150814871208e-43 |
p.valor <- boxM(data = df_svm[, 1:12], grouping = df_svm$variedad)$p.value
Test <- boxM(data = df_svm[, 1:12], grouping = df_svm$variedad)$method
# Imprimo resultados de test
kable(cbind(Test, p.valor))
| Test | p.valor |
|---|---|
| Box’s M-test for Homogeneity of Covariance Matrices | 0 |
matriz_de_distancias <- vegan::betadisper(dist(df_svm[, 1:12], method='euclidean'), df_svm$variedad, type = c("median","centroid"), bias.adjust = T,sqrt.dist = FALSE, add = FALSE)
## Registered S3 methods overwritten by 'vegan':
## method from
## rev.hclust dendextend
## plot.rda klaR
## predict.rda klaR
## print.rda klaR
test_levene <- anova(matriz_de_distancias)
p.valor <- test_levene$`Pr(>F)`[1]
TukeyHSD(matriz_de_distancias)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = distances ~ group, data = df)
##
## $group
## diff lwr upr p adj
## 2-1 -12.98251 -15.28052 -10.68451 0
#plot(matriz_de_distancias)
Test <- 'Levene'
# Imprimo resultado del test
kable(cbind(Test,p.valor))
| Test | p.valor |
|---|---|
| Levene | 9.94794504885316e-28 |
HOTELLING <- HotellingsT2Test(as.matrix(df_svm[,-c(13)]) ~ variedad, data =df_svm)
Test <- HOTELLING$method
p.valor <- HOTELLING$p.value[1]
# Imprimo resultados en una tabla
kable(cbind(Test, p.valor ))
| Test | p.valor |
|---|---|
| Hotelling’s two sample T2-test | 0 |
# se utiliza el paquete npmv no paramétrico para comparar vector de medias. (Nonparametric Inference for Multivariate Data: R Package npmv, January 2017, Volume 76, Issue 4. doi: 10.18637/jss.v076.i04, https://www.jstatsoft.org/article/view/v076i04)
# => Agregar todas las variables
#noparam <- nonpartest(pH | calidad | densidad | alcohol | sulfatos | cloruros ~ variedad, data = df_svm, permreps = 1000, plots=F)
noparam <- nonpartest( acidez_fija | acidez_volatil | acido_citrico | azucar_residual | cloruros | anhidrido_sulfuroso_libre | anhidrido_sulfuroso_total | densidad | pH | sulfatos | alcohol | calidad ~ variedad, data = df_svm, permreps = 1000, plots=F)
p.valor <- noparam$results$`P-value`[1]
Test <- 'No paramétrico multivariado (npmv)'
# Imprimo el resultado
kable(cbind(Test, p.valor ))
| Test | p.valor |
|---|---|
| No paramétrico multivariado (npmv) | 0 |
set.seed(661) # para asegurar reproducibilidad
dt = sort(sample(nrow(df_svm), nrow(df_svm)*.8))
datos_train_svm<-df_svm[dt,]
datos_test_svm<-df_svm[-dt,]
# Se realiza el escalado/estandarización con -> (x - mean(x)) / sd(x)
# Calculo media y sd de subconjunto de entrenamiento (train), y con esos datos hago el escalado del test. La idea de escalar el conjunto de test (prueba) utilizando datos solamente de train es para evitar el data leakeage.
# Hago escalado a mano del test set con media del training, y sd del training
for (k in 1:12){
datos_test_svm[,k]=round(print(datos_test_svm[,..k]-colMeans(datos_test_svm[,..k]))/sd(unlist(datos_train_svm[,..k])),2)
}
## acidez_fija
## 1: -0.53625
## 2: -0.33625
## 3: -1.13625
## 4: -1.13625
## 5: -2.73625
## ---
## 396: -0.43625
## 397: -0.03625
## 398: 3.16375
## 399: -1.63625
## 400: -0.83625
## acidez_volatil
## 1: -0.1621625
## 2: -0.1821625
## 3: 0.2278375
## 4: -0.0521625
## 5: -0.0421625
## ---
## 396: -0.0721625
## 397: 0.3128375
## 398: -0.1821625
## 399: 0.1678375
## 400: 0.0178375
## acido_citrico
## 1: -0.12185
## 2: -0.04185
## 3: -0.11185
## 4: 0.02815
## 5: -0.06185
## ---
## 396: -0.08185
## 397: -0.21185
## 398: 0.17815
## 399: -0.31185
## 400: -0.00185
## azucar_residual
## 1: 9.123625
## 2: -3.376375
## 3: 4.123625
## 4: -1.976375
## 5: 3.223625
## ---
## 396: -2.276375
## 397: -1.976375
## 398: -1.776375
## 399: -2.376375
## 400: -2.576375
## cloruros
## 1: -0.014905
## 2: -0.025905
## 3: -0.020905
## 4: -0.015905
## 5: -0.035905
## ---
## 396: -0.000905
## 397: 0.017095
## 398: 0.021095
## 399: 0.012095
## 400: 0.012095
## anhidrido_sulfuroso_libre
## 1: 21.765
## 2: 1.765
## 3: 55.765
## 4: 15.265
## 5: -1.235
## ---
## 396: 9.765
## 397: -16.235
## 398: -14.235
## 399: -19.235
## 400: -4.235
## anhidrido_sulfuroso_total
## 1: 102.85625
## 2: 4.85625
## 3: 129.85625
## 4: 115.85625
## 5: 21.85625
## ---
## 396: -24.14375
## 397: -68.14375
## 398: -62.14375
## 399: -85.14375
## 400: -43.14375
## densidad
## 1: 0.00306995
## 2: -0.00288005
## 3: 0.00046995
## 4: -0.00244005
## 5: -0.00317005
## ---
## 396: 0.00029995
## 397: 0.00201995
## 398: 0.00161995
## 399: -0.00078005
## 400: 0.00109995
## pH
## 1: -0.1177
## 2: -0.0477
## 3: -0.0677
## 4: -0.0877
## 5: 0.1523
## ---
## 396: 0.1923
## 397: 0.1523
## 398: -0.0177
## 399: 0.3523
## 400: 0.2323
## sulfatos
## 1: -0.033625
## 2: -0.233625
## 3: -0.123625
## 4: 0.046375
## 5: -0.163625
## ---
## 396: 0.056375
## 397: -0.073625
## 398: 0.116375
## 399: -0.023625
## 400: -0.013625
## alcohol
## 1: -1.4245417
## 2: -0.6245417
## 3: -1.1245417
## 4: -0.4245417
## 5: 0.8754583
## ---
## 396: -0.3245417
## 397: -0.7245417
## 398: 1.2754583
## 399: 1.0754583
## 400: -0.9245417
## calidad
## 1: -0.7375
## 2: 0.2625
## 3: -0.7375
## 4: -0.7375
## 5: 0.2625
## ---
## 396: -0.7375
## 397: -0.7375
## 398: 0.2625
## 399: 0.2625
## 400: 0.2625
# Hago automáticamente con la función scale, el escalado de training
datos_train_svm = as.data.frame(scale(datos_train_svm[,1:12]))
datos_train_svm$variedad <- df_svm[dt,]$variedad
# y las pongo en = orden que testset
# creo una única df con todos los datos escalados (que usaré luego para clustering)
datos_escalados <- as.data.frame(scale(df_svm[,1:12]))
datos_escalados$variedad <- df_svm$variedad
# Defino modelo SVM => Tengo que agregar todas las columnas
sample_svn_tr = datos_train_svm[,1:13]
sample_svn_te = datos_test_svm[,1:13]
set.seed(661)
task = makeClassifTask(data = sample_svn_tr, target = "variedad")
lrn_svm1 = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "linear", cost=2))
mod_svm1 = mlr::train(lrn_svm1, task)
# Predicción TEST
pred_svm1= predict(mod_svm1, newdata = sample_svn_te)
## Warning in predict.WrappedModel(mod_svm1, newdata = sample_svn_te): Provided
## data for prediction is not a pure data.frame but from class data.table, hence it
## will be converted.
acc_svm1 <- round(measureACC(pred_svm1$data$truth, pred_svm1$data$response),3)
AUC_svm1_te <- round(measureAUC(as.data.frame(pred_svm1)$prob.2,as.data.frame(pred_svm1)$truth,1,2),3)
# ···························································································
# Predicción TRAIN (naive)
pred_svm_1 = predict(mod_svm1, newdata = sample_svn_tr) # para comparar con set de test
acc_svm_1 <- round(measureACC(as.data.frame(pred_svm_1)$truth, as.data.frame(pred_svm_1)$response),3)
AUC_svm1_tr <- round(measureAUC(as.data.frame(pred_svm_1)$prob.2,as.data.frame(pred_svm_1)$truth, 1,2),3)
#print(sample_svn_tr)
# ···························································································
# Cambio el threshold [esto lo hago para train y test]
acc=NULL
acc2=NULL
threshold = seq(0.1,0.95,0.01)
for (i in 1:length(threshold)) {
pred = setThreshold(pred_svm1, threshold = threshold[i])
acc[i] = measureACC(as.data.frame(pred)$truth, as.data.frame(pred)$response)}
for (i in 1:length(threshold)) {
pred2 = setThreshold(pred_svm_1, threshold = threshold[i])
acc2[i] = measureACC(as.data.frame(pred2)$truth, as.data.frame(pred2)$response)}
par(mfrow = c(2,5))
new_df1 <- as.data.frame(cbind(threshold,acc))
new_df1 <- new_df1%>%mutate(sub_data='test')
new_df2 <- as.data.frame(cbind(threshold,acc2))
colnames(new_df2) <- c('threshold','acc')
new_df2 <- new_df2%>%mutate(sub_data='train')
new_df <- as.data.frame(rbind(new_df1,new_df2))
#new_df1[which.max(new_df1$acc),"threshold"]
#new_df2[which.max(new_df2$acc),"threshold"]
# ···························································································
# Gráfico de cómo varía la métrica de performance accuracy, de acuerdo al umbral elegido
ggplot(new_df, aes(x=threshold, y=acc)) + geom_line(aes(color = sub_data,linetype=sub_data)) +
theme +labs(x='Umbral', y='Métrica de performance (accuracy)',
title= 'Evaluación del modelo de Máquinas de soporte vectorial SVM') +
scale_color_manual(values = c("red", "darkred"),labels=c('prueba','entrenamiento')) +
scale_linetype_manual(values=c(1,2), labels=c('prueba','entrenamiento')) +
labs(color='Conjunto de\n evaluación',linetype='Conjunto de\n evaluación')
#Representamos la matriz de confusion para poder hacer un analisis de las diferentes metricas sobre la prediccion
#install.packages("caret")
library(caret)
##
## Attaching package: 'caret'
## The following objects are masked from 'package:DescTools':
##
## MAE, RMSE
## The following object is masked from 'package:mlr':
##
## train
## The following object is masked from 'package:purrr':
##
## lift
confusionMatrix(as.data.frame(pred_svm1)$truth, as.data.frame(pred_svm1)$response)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 1 2
## 1 199 0
## 2 2 199
##
## Accuracy : 0.995
## 95% CI : (0.9821, 0.9994)
## No Information Rate : 0.5025
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.99
##
## Mcnemar's Test P-Value : 0.4795
##
## Sensitivity : 0.9900
## Specificity : 1.0000
## Pos Pred Value : 1.0000
## Neg Pred Value : 0.9900
## Prevalence : 0.5025
## Detection Rate : 0.4975
## Detection Prevalence : 0.4975
## Balanced Accuracy : 0.9950
##
## 'Positive' Class : 1
##
df_svm_threshold = generateThreshVsPerfData(list(svm_te = pred_svm1, svm_tr = pred_svm_1),
measures = list(fpr, tpr, mmce))
plotROCCurves(df_svm_threshold) + theme +
labs(title='Curva ROC del modelo de Máquinas de soporte vectorial SVM kernel lineal',
x='Tasa de falsos positivos (FPR)', y='Tasa de positivos verdaderos (TPR)',
color='Conjunto de\n evaluación') +
scale_color_manual(values = c("red", "darkred"), labels=c('prueba','entrenamiento')) +
geom_label(label='AUC test', x=0.35, y=0.75, label.size = 0.3, size=4,
color = "red",fill="white") +
geom_label(label='AUC train',x=0.07, y=0.97, label.size = 0.3, size=4,
color = "darkred",fill="white")
# ················ Métricas del modelo de SVM ················
Métrica <- c('valor','datos')
Accuracy <- c(acc_svm1,'prueba')
Accuracy. <- c(acc_svm_1,'entrenamiento')
AUC_ROC <- c(AUC_svm1_te,'prueba')
AUC_ROC. <- c(AUC_svm1_tr,'entrenamiento')
# Imprimo resultados
kable(rbind(Métrica, Accuracy, Accuracy., AUC_ROC, AUC_ROC.))
| Métrica | valor | datos |
| Accuracy | 0.995 | prueba |
| Accuracy. | 0.992 | entrenamiento |
| AUC_ROC | 0.998 | prueba |
| AUC_ROC. | 0.996 | entrenamiento |
# Defino modelo SVM radial
sample_svn_tr_radial = datos_train_svm[,1:13]
sample_svn_te_radial = datos_test_svm[,1:13]
set.seed(661)
task = makeClassifTask(data = sample_svn_tr_radial, target = "variedad")
lrn_svm3 = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "radial", cost=2))
mod_svm3 = mlr::train(lrn_svm3, task)
# Predicción TEST
pred_svm3= predict(mod_svm3, newdata = sample_svn_te_radial)
## Warning in predict.WrappedModel(mod_svm3, newdata = sample_svn_te_radial):
## Provided data for prediction is not a pure data.frame but from class data.table,
## hence it will be converted.
acc_svm3 <- round(measureACC(pred_svm3$data$truth, pred_svm3$data$response),3)
AUC_svm3_te <- round(measureAUC(as.data.frame(pred_svm3)$prob.2,as.data.frame(pred_svm3)$truth,1,2),3)
# ···························································································
# Predicción TRAIN (naive)
pred_svm_3 = predict(mod_svm3, newdata = sample_svn_tr_radial) # por si quiero ver naive sobre training
acc_svm_3 <- round(measureACC(as.data.frame(pred_svm_3)$truth, as.data.frame(pred_svm_3)$response),3)
AUC_svm3_tr <- round(measureAUC(as.data.frame(pred_svm_3)$prob.2,as.data.frame(pred_svm_3)$truth, 1,2),3)
# ···························································································
# Cambio el threshold [esto lo hago para train y test]
acc=NULL
acc2=NULL
threshold = seq(0.1,0.95,0.01)
for (i in 1:length(threshold)) {
pred = setThreshold(pred_svm3, threshold = threshold[i])
acc[i] = measureACC(as.data.frame(pred)$truth, as.data.frame(pred)$response)}
for (i in 1:length(threshold)) {
pred2 = setThreshold(pred_svm_3, threshold = threshold[i])
acc2[i] = measureACC(as.data.frame(pred2)$truth, as.data.frame(pred2)$response)}
par(mfrow = c(2,5))
new_df1 <- as.data.frame(cbind(threshold,acc))
new_df1 <- new_df1%>%mutate(sub_data='test')
new_df2 <- as.data.frame(cbind(threshold,acc2))
colnames(new_df2) <- c('threshold','acc')
new_df2 <- new_df2%>%mutate(sub_data='train')
new_df <- as.data.frame(rbind(new_df1,new_df2))
#new_df1[which.max(new_df1$acc),"threshold"]
#new_df2[which.max(new_df2$acc),"threshold"]
# ···························································································
# Gráfico de cómo varía la métrica de performance accuracy, de acuerdo al umbral elegido
ggplot(new_df, aes(x=threshold, y=acc)) + geom_line(aes(color = sub_data,linetype=sub_data)) +
theme +labs(x='Umbral', y='Métrica de performance (accuracy)',
title= 'Evaluación del modelo de Máquinas de soporte vectorial SVM') +
scale_color_manual(values = c("red", "darkred"),labels=c('prueba','entrenamiento')) +
scale_linetype_manual(values=c(1,2), labels=c('prueba','entrenamiento')) +
labs(color='Conjunto de\n evaluación',linetype='Conjunto de\n evaluación')
# Para independizarnos de la elección del umbral, grafico curvas ROC para las predicciones del modelo SVM con los datos de TEST y TRAIN
df_svm_threshold_radial = generateThreshVsPerfData(list(svm_te = pred_svm3, svm_tr = pred_svm_3),
measures = list(fpr, tpr, mmce))
plotROCCurves(df_svm_threshold_radial)
confusionMatrix(as.data.frame(pred_svm3)$truth, as.data.frame(pred_svm3)$response)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 1 2
## 1 199 0
## 2 2 199
##
## Accuracy : 0.995
## 95% CI : (0.9821, 0.9994)
## No Information Rate : 0.5025
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.99
##
## Mcnemar's Test P-Value : 0.4795
##
## Sensitivity : 0.9900
## Specificity : 1.0000
## Pos Pred Value : 1.0000
## Neg Pred Value : 0.9900
## Prevalence : 0.5025
## Detection Rate : 0.4975
## Detection Prevalence : 0.4975
## Balanced Accuracy : 0.9950
##
## 'Positive' Class : 1
##
# ················ Métricas del modelo de SVM ················
Métrica <- c('valor','datos')
Accuracy <- c(acc_svm3,'prueba')
Accuracy. <- c(acc_svm_3,'entrenamiento')
AUC_ROC <- c(AUC_svm3_te,'prueba')
AUC_ROC. <- c(AUC_svm3_tr,'entrenamiento')
# Imprimo resultados
kable(rbind(Métrica, Accuracy, Accuracy., AUC_ROC, AUC_ROC.))
| Métrica | valor | datos |
| Accuracy | 0.995 | prueba |
| Accuracy. | 0.995 | entrenamiento |
| AUC_ROC | 0.996 | prueba |
| AUC_ROC. | 0.999 | entrenamiento |
SVM_metrics1 <- calculateROCMeasures(pred_svm1)
SVM_metrics3 <- calculateROCMeasures(pred_svm3)
# Ingenua para ver cómo le va
pred_todos=NULL
pred_todos_svm1 <- as.data.frame(predict(mod_svm1, newdata = datos_escalados[-13]))
pred_todos_svm3 <- as.data.frame(predict(mod_svm3, newdata = datos_escalados[-13]))
# ······················ PCA datos originales ··················································
ggbiplot(df_svm.pc, obs.scale=0.1 ,var.scale=1,alpha=0.5,groups=factor(df_pca$variedad)) +
scale_color_manual(name="variedad", values=c('royalblue2','#ff7474ff'),
labels=c("blanco","tinto")) +
theme + labs(
title= 'Representación de las etiquetas reales\nen las componentes principales 1 y 2')+
theme(legend.position=c(.865,.15))
# ······················ PCA + predicciones SVM lineal (modelo 1) ·······························
ggbiplot(df_svm.pc, obs.scale=0.1 ,var.scale=1,alpha=.5,groups=factor(pred_todos_svm1$response)) +
scale_color_manual(name="Predicción", values=c('royalblue2','#ff7474ff'),
labels=c("blanco","tinto")) +
theme + labs(title= 'Representación de las predicciones ingenuas del modelo SVM (l)\nen las componentes principales 1 y 2') +
theme(legend.position=c(.865,.15))
# ······················ PCA + predicciones SVM radial (modelo 3) ·······························
ggbiplot(df_svm.pc, obs.scale=0.1 ,var.scale=1,alpha=.5,groups=factor(pred_todos_svm3$response)) +
scale_color_manual(name="Predicción", values=c('royalblue2','#ff7474ff'),
labels=c("blanco","tinto")) +
theme + labs(title= 'Representación de las predicciones ingenuas del modelo SVM (r)\nen las componentes principales 1 y 2') +
theme(legend.position=c(.865,.15))
# ······················ curvas ROC para todos los modelos ········································
df_todos = generateThreshVsPerfData(list(svm1 = pred_svm1,
svm3 = pred_svm3), measures = list(fpr, tpr, mmce))
plotROCCurves(df_todos) + theme +
labs(title='Curvas ROC de modelos de clasificación supervisada (datos de prueba)',
x='Tasa de falsos positivos (FPR)', y='Tasa de positivos verdaderos (TPR)',
color=' Modelo en\n evaluación') +
scale_color_manual(values = c("red", "black"),
labels=c('SVM (l)','SVM (r)'))+
theme(legend.position=c(0.915,0.25))
Conclusiones Podemos determinar a partir del análisis realizado, de que para este caso, SVM está realizando una separación muy precisa de los datos llegando a una precisión por encima del 99% tanto en el conjunto de prueba, como de validación
Como mencionamos anteriormente los supuestos no se cumplen por lo que no podemos asegurar que los resultados sean totalmente confiables
Por otro lado, los porcentajes de comparacion entre prueba y entrenamiento son sospechosamente altos, por lo que podriamos pensar que se esta generando un sobreajuste a pesar de separar el conjunto de datos en entretamiento y prueba
Sería adecuado poder hacer este análisis con un conjunto de datos mucho mayor, para poder tener una mayor confiabilidad en el modelo
Casualmente dan iguales o mejores resultados en el conjunto de prueba que en el de entrenamiento, considero que esto se da por la alta precisión del modelo y la baja cantidad de datos para probar
df_clustering = sample_data
cantidad_clusters=2
set.seed(661)
par(mfcol = c(1,1))
data_c_diag = df_clustering[-c(13)] #
sample_cluster1 <- data_c_diag[sample(1:nrow(data_c_diag), 1000,replace=FALSE),]
sample_cluster <- as.data.frame(scale(sample_cluster1[,1:12]))
sample_cluster$variedad <- sample_cluster1$variedad
# Imprimo cuántos valores hay en cada categoría de la variable diagnóstico en mi muestra de n=100
kable(table(sample_cluster$variedad))
| Var1 | Freq |
|---|---|
| 1 | 503 |
| 2 | 497 |
# Escalo los datos y hago PCA
df_sample_cluster = sample_cluster[-13]
df_sample_cluster_pca <-df_sample_cluster %>% dplyr::select(is.numeric)
df_clustering.pc2 = prcomp(df_sample_cluster_pca,scale = TRUE)
# Matriz de distancias euclídeas
mat_dist <- dist(x = df_sample_cluster, method = "euclidean")
# Dendrogramas (según el tipo de segmentación jerárquica aplicada)
hc_complete <- hclust(d = mat_dist, method = "complete")
hc_average <- hclust(d = mat_dist, method = "average")
hc_single <- hclust(d = mat_dist, method = "single")
hc_ward <- hclust(d = mat_dist, method = "ward.D2")
#Calculo del coeficiente de correlacion cofenetico
completo <- round(cor(x = mat_dist, cophenetic(hc_complete)),3)
promedio <- round(cor(x = mat_dist, cophenetic(hc_average)),3)
simple <- round(cor(x = mat_dist, cophenetic(hc_single)),3)
ward <- round(cor(x = mat_dist, cophenetic(hc_ward)),3)
valores_coef <- cbind(completo,promedio,simple,ward)
# Imprimo valores de coeficiente cofenético
kable(valores_coef)
| completo | promedio | simple | ward |
|---|---|---|---|
| 0.587 | 0.772 | 0.649 | 0.431 |
# Armo clusters
jer_ward<-cutree(hc_ward,k=cantidad_clusters)
jer_average<-cutree(hc_average,k=cantidad_clusters)
jer_complete<-cutree(hc_complete,k=cantidad_clusters)
jer_single<-cutree(hc_single,k=cantidad_clusters)
# Agrego cluster a dataframe
sample_cluster$jer_ward=jer_ward
sample_cluster$jer_average=jer_average
sample_cluster$jer_complete=jer_complete
sample_cluster$jer_single=jer_single
mar = c(5.1, 4.1, 4.1, 2.1)
pch=c('royalblue2','#ff7474ff')
cols=alpha(pch[sample_cluster$variedad[order.dendrogram(as.dendrogram(hc_ward))]],0.7)
dend_ward <- color_branches(as.dendrogram(hc_ward), k = 2)
dend_ward <- set(dend_ward, "labels_cex", 0.1)
grafico1 <- dend_ward %>% set("leaves_pch",19)%>%
set("leaves_cex", .9) %>% set("leaves_col", cols) %>%
plot(main = "Dendrograma jerárquico", ylab='Distancia',cex.lab=1, cex.axis=.6)+
mtext(side = 3, line = 0.5, at = 1, adj = -1.8, 'Distancia Ward')+
mtext(side = 1, line = 0.5, at = 1, adj = -4.1, 'Vinos')
legend(5,28, title='Variedad',
legend = c("blanco" , "tinto"),
col = c('royalblue2','#ff7474ff') ,
pch = c(19,19), bty = "n", pt.cex = 1.5, cex = 0.8 ,
text.col = "black", horiz = FALSE, inset = c(0, 0.1))
cols_a=alpha(pch[sample_cluster$variedad[order.dendrogram(as.dendrogram(hc_average))]],0.7)
dend_average <- color_branches(as.dendrogram(hc_average), k = 2)
dend_average <- set(dend_average, "labels_cex", 0.1)
grafico2 <- dend_average %>% set("leaves_pch",19) %>%
set("leaves_cex", .8) %>% set("leaves_col", cols_a) %>%
plot(main = "Dendrograma jerárquico", ylab='Distancia',cex.lab=1, cex.axis=.6)+
mtext(side = 3, line = 0.5, at = 1, adj = -1.3, 'Distancia Promedio')+
mtext(side = 1, line = 0.5, at = 1, adj = -4.1, 'Vinos')
legend(76,8, title='Variedad',
legend = c("blanco" , "tinto"),
col = c('royalblue2','#ff7474ff') ,
pch = c(19,19), bty = "n", pt.cex = 1.5, cex = 0.8 ,
text.col = "black", horiz = FALSE, inset = c(0, 0.1))
cols_c=alpha(pch[sample_cluster$variedad[order.dendrogram(as.dendrogram(hc_complete))]],0.7)
dend_complete <- color_branches(as.dendrogram(hc_complete), k = 2)
dend_complete <- set(dend_complete, "labels_cex", 0.1)
grafico3 <- dend_complete %>% set("leaves_pch",19) %>% # node point type
set("leaves_cex", .8) %>% # node point size
set("leaves_col", cols_c) %>% #node point color
plot(main = "Dendrograma jerárquico", ylab='Distancia',cex.lab=1, cex.axis=.6)+
mtext(side = 3, line = 0.5, at = 1, adj = -1.3, 'Distancia Completa')+
mtext(side = 1, line = 0.5, at = 1, adj = -4.1, 'Vinos')
legend(85,16, title='Diagnóstico',
legend = c("blanco" , "tinto"),
col = c('royalblue2','#ff7474ff') ,
pch = c(19,19), bty = "n", pt.cex = 1.5, cex = 0.8 ,
text.col = "black", horiz = FALSE, inset = c(0, 0.1))
cols_s=alpha(pch[sample_cluster$variedad[order.dendrogram(as.dendrogram(hc_single))]],0.7)
dend_single <- color_branches(as.dendrogram(hc_single), k = 2)
dend_single <- set(dend_single, "labels_cex", 0.1)
grafico4 <- dend_single %>% set("leaves_pch",19) %>% # node point type
set("leaves_cex", .8) %>% # node point size
set("leaves_col", cols_s) %>% #node point color
plot(main = "Dendrograma jerárquico", ylab='Distancia',cex.lab=1, cex.axis=.7)+
mtext(side = 3, line = 0.5, at = 1, adj = -1.7, 'Distancia Simple')+
mtext(side = 1, line = 0.5, at = 1, adj = -4.1, 'Vinos')
legend(80,7, title='Variedad',
legend = c("blanco" , "tinto"),
col = c('royalblue2','#ff7474ff') ,
pch = c(19,19), bty = "n", pt.cex = 1.5, cex = 0.8 ,
text.col = "black", horiz = FALSE, inset = c(0, 0.1))
# ·····················································
# cuántos pacientes de cada diagnóstico están en cada cluster:
ward_cluster1 <- sample_cluster %>% filter (jer_ward == '1')
cluster1 <- table(ward_cluster1$variedad)
Ward_cluster.1 <- round(prop.table(cluster1)*100,2)
# ·····················································
ward_cluster2 <- sample_cluster %>% filter (jer_ward == '2')
cluster2 <- table(ward_cluster2$variedad)
Ward_cluster.2 <- round(prop.table(cluster2)*100,2)
kable(cbind(rbind(cluster1,cluster2),rbind(Ward_cluster.1,Ward_cluster.2)))
| 1 | 2 | 1 | 2 | |
|---|---|---|---|---|
| cluster1 | 32 | 490 | 6.13 | 93.87 |
| cluster2 | 471 | 7 | 98.54 | 1.46 |
# ·····················································
promedio_cluster1 <- sample_cluster %>% filter (jer_average == '1')
cluster1 <- table(promedio_cluster1$variedad)
promedio_cluster.1 <- round(prop.table(cluster1)*100,2)
# ·····················································
promedio_cluster2 <- sample_cluster %>% filter (jer_average == '2')
cluster2 <- table(promedio_cluster2$variedad)
promedio_cluster.2 <- round(prop.table(cluster2)*100,2)
kable(cbind(rbind(cluster1,cluster2),rbind(promedio_cluster.1,promedio_cluster.2)))
| 1 | 2 | 1 | 2 | |
|---|---|---|---|---|
| cluster1 | 502 | 497 | 50.25 | 49.75 |
| cluster2 | 1 | 0 | 100.00 | 0.00 |
# ·····················································
completo_cluster1 <- sample_cluster %>% filter (jer_complete == '1')
cluster1 <- table(completo_cluster1$variedad)
completo_cluster.1 <- round(prop.table(cluster1)*100,2)
# ·····················································
completo_cluster2 <- sample_cluster %>% filter (jer_complete == '2')
cluster2 <- table(completo_cluster2$variedad)
completo_cluster.2 <- round(prop.table(cluster2)*100,2)
kable(cbind(rbind(cluster1,cluster2),rbind(completo_cluster.1,completo_cluster.2)))
| 1 | 2 | 1 | 2 | |
|---|---|---|---|---|
| cluster1 | 502 | 497 | 50.25 | 49.75 |
| cluster2 | 1 | 0 | 100.00 | 0.00 |
# ·····················································
simple_cluster1 <- sample_cluster %>% filter (jer_single == '1')
cluster1 <- table(simple_cluster1$variedad)
simple_cluster.1 <- round(prop.table(cluster1)*100,2)
# ·····················································
simple_cluster2 <- sample_cluster %>% filter (jer_single == '2')
cluster2 <- table(simple_cluster2$variedad)
simple_cluster.2 <- round(prop.table(cluster2)*100,2)
kable(cbind(rbind(cluster1,cluster2),rbind(simple_cluster.1,simple_cluster.2)))
| 1 | 2 | 1 | 2 | |
|---|---|---|---|---|
| cluster1 | 502 | 497 | 50.25 | 49.75 |
| cluster2 | 1 | 0 | 100.00 | 0.00 |
Conclusiones Como resultado del análisis de clasificación realizada utilizando diferentes tipos de segmentación jerárquica, determinamos que a pesar de que el cálculo del coeficiente de correlación cofenético nos arrojo con mayor valor al tipo promedio, ejecutando el denograma de cada uno vemos que claramente el que mejor realiza la separación de clusteres es el tipo ward, arrojando resultados muchísimos mejores que el resto
1 2 1 2 cluster1 32 490 6.13 93.87 cluster2 471 7 98.54 1.46
Limito el numero de clusters por el costo computacional
#datos_para_kmeans = sample_data[1:12]
df_kmeans = sample_data
datos_para_kmeans <- sample_data %>% dplyr::select(is.numeric)
#analisis de la cantidad de clusters. Este primer bloque es solo para definir funciones.
#se define una funcion para calcular metricas que orientan sobre el numero de clusters a elegir para el problema.
metrica = function(datA_esc,kmax,f) {
sil = array()
sse = array()
datA_dist= dist(datA_esc,method = "euclidean", diag = FALSE, upper = FALSE, p = 2)
for ( i in 2:kmax) { if (strcmp(f,"kmeans")==TRUE) { #centroide: tipico kmeans
CL = kmeans(datA_esc,centers=i,nstart=50,iter.max = kmax)
sse[i] = CL$tot.withinss
CL_sil = silhouette(CL$cluster, datA_dist)
sil[i] = summary(CL_sil)$avg.width}
if (strcmp(f,"pam")==TRUE){ #medoide: ojo porque este metodo tarda muchisimo
CL = pam(x=datA_esc, k=i, diss = F, metric = "euclidean")
sse[i] = CL$objective[1]
sil[i] = CL$silinfo$avg.width}}
sse
sil
return(data.frame(sse,sil))}
#en este bloque se estudia cuantos clusters convendría generar segun indicadores tipicos -> por ejemplo el "Silhouette"
kmax = 10
m1 = metrica(scale(datos_para_kmeans),kmax,"kmeans") #tipica con estimadores de la normal
m1 <- m1[complete.cases(m1),]
m1$kcluster <- seq(2,kmax,1)
m1 <- m1%>%dplyr::select(3,1,2)
m1_sse <- m1%>%dplyr::select(-3)%>%mutate(metric='SSE')
colnames(m1_sse) <- c('kcluster','value','metric')
m1_sil <- m1%>%dplyr::select(-2)%>%mutate(metric='SIL')
colnames(m1_sil) <- c('kcluster','value','metric')
m1 <- rbind(m1_sse,m1_sil)
# Grafico de métricas SIL y SSE
ggplot(m1, aes(kcluster, value, linetype=metric)) + geom_line(col='red') +
facet_wrap(~metric, ncol=1, scales='free')+theme+geom_point(col='red', size=2, fill='pink', shape=21)+
labs(title='Determinación de número de clusters',
x='k Número de clusters', y='Valor', linetype='Métrica')+
scale_x_continuous(breaks = seq(1, kmax, by = 1))+
scale_linetype_manual(values=c(1,2))
set.seed(661)
cantidad_clusters=2
CL = kmeans(scale(datos_para_kmeans),cantidad_clusters)
df_kmeans$kmeans = CL$cluster
# Grafico scatterplot original + cluster con k=2
par(mfrow=c(1,2))
# -------------------------------------------------------------------------------------------
col1 <- c('royalblue2','#ff7474ff')
col1 <- col1[as.numeric(df_kmeans$variedad)]
scatterplot3d(df_kmeans$pH,df_kmeans$densidad,df_kmeans$calidad, color = alpha(col1,0.3), box=F,angle=45, pch = 19, grid = TRUE, tick.marks = FALSE, xlab = "ph", ylab = "densiidad", zlab = "calidad", main='Realidad')
#legend("topright", bty = "n", cex = .9, title = "Variedad", c("blanco", "tinto"), fill = c('royalblue2','#ff7474ff'))
# -------------------------------------------------------------------------------------------
colors <- c('orange','#a25da2a5')
colors <- colors[as.numeric(df_kmeans$kmeans)]
scatterplot3d(df_kmeans$pH,df_kmeans$densidad,df_kmeans$calidad, color = alpha(colors,0.3), box=F,angle=45, pch = 19, grid = TRUE, tick.marks = FALSE, xlab = "ph", ylab = "densiidad", zlab = "calidad", main='Clustering')
#legend("topright", bty = "n", cex = .9, title = "Grupo k-means", c("1", "2"), fill = c('orange','#a25da2a5'))
#conviene en un biplot ya que tengo las flechas de las variables originales
# GRAFICO ORIGINAL
ggbiplot(df_svm.pc, obs.scale=0.1 ,var.scale=1,alpha=0.5,
groups=factor(df_kmeans$variedad)) +
scale_color_manual(name="variedad",values=c('royalblue2','#ff7474ff'),
labels=c("blanco",'tinto')) +
theme + labs(title='Análisis de componentes principales') +
theme(legend.position=c(.85,.15)) +
labs(x='PC1 (51.5% explicado)', y= 'PC2 (22.5% explicado)')
# -------------------------------------------------------------------------------------------
# GRAFICO KMEANS
ggbiplot(df_svm.pc, obs.scale=0.1 ,var.scale=1, alpha=0.5,groups = as.factor(df_kmeans$kmeans) )+
scale_color_manual(name="Grupo k-means", values=c("orange",'#a25da2a5',"darkgreen"),
labels=c("1", "2","3")) + theme+
labs(title= 'Representación del clustering utilizando k-means')+theme(legend.position=c(.85,.15))
# ·····················································
# cuántos pacientes de cada diagnóstico están en cada cluster:
pacientes_cluster1 <- df_kmeans %>% filter (kmeans == '1')
cluster1 <- table(pacientes_cluster1$variedad)
porcentaje_cluster.1 <- round(prop.table(cluster1)*100,2)
pacientes_cluster2 <- df_kmeans %>% filter (kmeans == '2')
cluster2 <- table(pacientes_cluster2$variedad)
porcentaje_cluster.2 <- round(prop.table(cluster2)*100,2)
# ·····················································
# Imprimo resultados
kable(cbind(rbind(cluster1,cluster2),rbind(porcentaje_cluster.1,porcentaje_cluster.2)))
| 1 | 2 | 1 | 2 | |
|---|---|---|---|---|
| cluster1 | 980 | 13 | 98.69 | 1.31 |
| cluster2 | 20 | 987 | 1.99 | 98.01 |
kmeans1 <- df_kmeans %>% filter(kmeans==1) %>%dplyr::select(c(1:12))%>% colMeans()
kmeans2 <- df_kmeans %>% filter(kmeans==2) %>%dplyr::select(c(1:12))%>% colMeans()
blanco <- df_kmeans %>%filter(variedad=='1') %>%dplyr::select(c(1:12))%>% colMeans()
tinto <- df_kmeans %>%filter(variedad=='2') %>%dplyr::select(c(1:12))%>% colMeans()
# Imprimo resultados
kable(rbind(kmeans1,blanco,kmeans2,tinto))
| acidez_fija | acidez_volatil | acido_citrico | azucar_residual | cloruros | anhidrido_sulfuroso_libre | anhidrido_sulfuroso_total | densidad | pH | sulfatos | alcohol | calidad | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| kmeans1 | 6.857150 | 0.2776788 | 0.3430111 | 6.538872 | 0.0459849 | 36.41138 | 141.04632 | 0.9940956 | 3.182699 | 0.4889325 | 10.49300 | 5.905337 |
| blanco | 6.853950 | 0.2820350 | 0.3403700 | 6.490900 | 0.0462090 | 35.94800 | 140.12700 | 0.9940928 | 3.185840 | 0.4907700 | 10.48465 | 5.871000 |
| kmeans2 | 8.340119 | 0.5239623 | 0.2740516 | 2.529146 | 0.0876207 | 15.59533 | 46.08937 | 0.9967512 | 3.311668 | 0.6625025 | 10.44361 | 5.612711 |
| tinto | 8.353700 | 0.5213300 | 0.2762100 | 2.549050 | 0.0876880 | 15.91300 | 46.34400 | 0.9967726 | 3.309430 | 0.6618800 | 10.45162 | 5.645000 |
# ya sabíamos del qqplot que las variables no son normales univariadas (si se separan por diagnóstico). por TLC, no obstante, uno puede aproximar las univariadas y hacer un test acorde para ver diferencias de medias
blanco_N <- df_kmeans%>%filter(variedad=='1')
tinto_M <- df_kmeans%>%filter(variedad=='2')
pH <- z.test(blanco_N$pH, sigma.x=sd(blanco_N$pH), tinto_M$pH, sigma.y=sd(tinto_M$pH), conf.level=0.95)$p.value
densidad <- z.test(blanco_N$densidad, sigma.x=sd(blanco_N$densidad), tinto_M$densidad, sigma.y=sd(tinto_M$densidad), conf.level=0.95)$p.value
calidad <- z.test(blanco_N$calidad, sigma.x=sd(blanco_N$calidad), tinto_M$calidad, sigma.y=sd(tinto_M$calidad), conf.level=0.95)$p.value
blanco_1 <- df_kmeans%>%filter(kmeans=='1')
tinto_2 <- df_kmeans%>%filter(kmeans=='2')
pH_c <- z.test(blanco_1$pH, sigma.x=sd(blanco_1$pH), tinto_2$pH, sigma.y=sd(tinto_2$pH), conf.level=0.95)$p.value
densidad_c <- z.test(blanco_1$densidad, sigma.x=sd(blanco_1$densidad), tinto_2$densidad, sigma.y=sd(tinto_2$densidad), conf.level=0.95)$p.value
calidad_c <- z.test(blanco_1$calidad, sigma.x=sd(blanco_1$calidad), tinto_2$calidad, sigma.y=sd(tinto_2$calidad), conf.level=0.95)$p.value
variable <- c('variedad','cluster')
# Imprimo resultados
kable(rbind(variable,cbind(rbind(pH, densidad, calidad),rbind(pH_c, densidad_c, calidad_c))))
| variable | variedad | cluster |
| pH | 7.24729192241589e-73 | 2.36190019262579e-80 |
| densidad | 5.95484451106308e-126 | 1.94113514456475e-122 |
| calidad | 5.17424960452661e-09 | 2.88274220391404e-14 |
#datos_para_kmeans = sample_data[1:12]
df_kmeans_4 = sample_data
datos_para_kmeans_4 <- sample_data %>% dplyr::select(is.numeric)
# ENTRENO K MEANS CON K=3
cantidad_clusters_4=4
set.seed(661)
CL_4 = kmeans(scale(datos_para_kmeans_4),cantidad_clusters_4)
df_kmeans_4$kmeans = CL_4$cluster
colors2 <- c('orange','#a25da2a5', 'darkgreen','lightblue')
colors2 <- colors2[as.numeric(df_kmeans_4$kmeans)]
# -------------------------------------------------------------------------------------------
# GRAFICO K=3
scatterplot3d(df_kmeans_4$pH,df_kmeans_4$densidad,df_kmeans_4$calidad, color = alpha(colors2,0.3), box=F,angle=25, pch = 19, grid = TRUE, tick.marks = FALSE, xlab = "pH", ylab = "densidad", zlab = "calidad", main='Clustering k=4', cex.lab=.8,scale.y=2, cex.symbols=1.3)
legend(x=4, y=5.5, bty = "n", cex = 1.1, title = "Grupo k-means", c("1", "2","3","4"), fill = c('orange','#a25da2a5', 'darkgreen','lightblue'))
# ·····················································
# cuántos pacientes de cada diagnóstico están en cada cluster:
pacientes_cluster1_4 <- df_kmeans_4 %>% filter (kmeans == '1')
cluster1_4 <- table(pacientes_cluster1_4$variedad)
porcentaje_cluster_4.1 <- round(prop.table(cluster1_4)*100,2)
pacientes_cluster2_4 <- df_kmeans_4 %>% filter (kmeans == '2')
cluster2_4 <- table(pacientes_cluster2_4$variedad)
porcentaje_cluster_4.2 <- round(prop.table(cluster2_4)*100,2)
pacientes_cluster3_4 <- df_kmeans_4 %>% filter (kmeans == '3')
cluster3_4 <- table(pacientes_cluster3_4$variedad)
porcentaje_cluster_4.3 <- round(prop.table(cluster3_4)*100,2)
pacientes_cluster4_4 <- df_kmeans_4 %>% filter (kmeans == '4')
cluster4_4 <- table(pacientes_cluster4_4$variedad)
porcentaje_cluster_4.4 <- round(prop.table(cluster4_4)*100,2)
# ·····················································
# Imprimo resultados
kable(cbind(rbind(cluster1_4,cluster2_4,cluster3_4,cluster4_4),rbind(porcentaje_cluster_4.1,porcentaje_cluster_4.2,porcentaje_cluster_4.3,porcentaje_cluster_4.4)))
| 1 | 2 | 1 | 2 | |
|---|---|---|---|---|
| cluster1_4 | 415 | 3 | 99.28 | 0.72 |
| cluster2_4 | 15 | 548 | 2.66 | 97.34 |
| cluster3_4 | 9 | 397 | 2.22 | 97.78 |
| cluster4_4 | 561 | 52 | 91.52 | 8.48 |
kmeans_41 <- df_kmeans_4 %>% filter(kmeans==1) %>%dplyr::select(c(1:12))%>% colMeans()
kmeans_42 <- df_kmeans_4 %>% filter(kmeans==2) %>%dplyr::select(c(1:12))%>% colMeans()
kmeans_43 <- df_kmeans_4 %>% filter(kmeans==3) %>%dplyr::select(c(1:12))%>% colMeans()
kmeans_44 <- df_kmeans_4 %>% filter(kmeans==4) %>%dplyr::select(c(1:12))%>% colMeans()
blanco <- df_kmeans_4 %>%filter(variedad=='1') %>%dplyr::select(c(1:12))%>% colMeans()
tinto <- df_kmeans_4 %>%filter(variedad=='2') %>%dplyr::select(c(1:12))%>% colMeans()
# Imprimo resultados
kable(rbind(kmeans_41,kmeans_42,kmeans_43,kmeans_44,blanco,tinto))
| acidez_fija | acidez_volatil | acido_citrico | azucar_residual | cloruros | anhidrido_sulfuroso_libre | anhidrido_sulfuroso_total | densidad | pH | sulfatos | alcohol | calidad | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| kmeans_41 | 6.927990 | 0.2887560 | 0.3575598 | 11.153469 | 0.0539545 | 45.71651 | 169.04426 | 0.9969552 | 3.157679 | 0.4869378 | 9.531738 | 5.543062 |
| kmeans_42 | 7.378686 | 0.6141297 | 0.1364121 | 2.399290 | 0.0808490 | 16.28686 | 49.15631 | 0.9963379 | 3.376927 | 0.5980462 | 10.196773 | 5.358792 |
| kmeans_43 | 9.816749 | 0.4133498 | 0.4635468 | 2.729926 | 0.0994852 | 14.54680 | 42.41626 | 0.9976136 | 3.208054 | 0.7437192 | 10.567406 | 5.884237 |
| kmeans_44 | 6.805791 | 0.2758401 | 0.3297227 | 3.129935 | 0.0414927 | 28.83524 | 115.68434 | 0.9921187 | 3.216444 | 0.5064600 | 11.290131 | 6.187602 |
| blanco | 6.853950 | 0.2820350 | 0.3403700 | 6.490900 | 0.0462090 | 35.94800 | 140.12700 | 0.9940928 | 3.185840 | 0.4907700 | 10.484650 | 5.871000 |
| tinto | 8.353700 | 0.5213300 | 0.2762100 | 2.549050 | 0.0876880 | 15.91300 | 46.34400 | 0.9967726 | 3.309430 | 0.6618800 | 10.451617 | 5.645000 |
Conclusiones Ejecutamos el modelo k means con el conjunto de datos, para esto en primer lugar realizamos el análisis correspondiente para elegir la cantidad de clusters a calcular. para esto calculamos las métricas Silhouette y SSE, obtenemos que el conjunto de datos seria bueno evaluarlo con 4, 5 clusters, sabiendo que nuestra variable target a dividir tan solo tiene dos categorías, continuamos la implementación con k means = 2 y k means = 4 para comparar resultados. Graficando un biplot con lo resuelto por K = 2 vemos que genera dos grupos muy marcados, donde no hay prácticamente un solapamiento entre ambos lo que sí sucede conjunto original, por el contrario con K = 4, si vemos que las fronteras entre los grupos comienzan a solaparse, por lo que con un análisis más profundo, la agrupación de varios clusters a través del análisis de sus medias podría brindar una imagen más acorde a la realidad. En este caso, observamos que k means = 1 y k means = 4 tienen medias muy similares, que podríamos asociarlas a la distribución de vino blanco.